Momentum distribution of the insulating phases of the extended Bose-Hubbard model 
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We develop two methods to calculate the momentum distribution of the insulating (Mott and charge-density- 
wave) phases of the extended Bose-Hubbard model with on-site and nearest-neighbor boson-boson repulsions 
on d-dimensional hypercubic lattices. First we construct the random phase approximation result, which corre- 
sponds to the exact solution for the infinite-dimensional limit. Then we perform a power-series expansion in the 
hopping t via strong-coupling perturbation theory, to evaluate the momentum distribution in two and three di- 
mensions; we also use the strong-coupling theory to verify the random phase approximation solution in infinite 
dimensions. Finally, we briefly discuss possible implications of our results in the context of ultracold dipolar 
Bose gases with dipole-dipole interactions loaded into optical lattices. 
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I. INTRODUCTION 

Ultracold atomic gases loaded into optical lattices have 
been proven to be ideal systems for studying Hubbard-type 
Hamiltonians JT|], the most successful of which has been the 
Bose-Hubbard (BH) model. This model has three terms J2]: 
a kinetic energy term which allows for the tunneling of the 
bosons between nearest-neighbor lattice sites, a potential en- 
ergy term which is given by the repulsion between bosons that 
occupy the same lattice site, and a chemical potential term 
which fixes the number of bosons. The phase diagram of this 
model has been known for a long time Q |3|, U, |5|, Is |7|] . The 
competition between the kinetic and potential energy terms 
leads to two phases: a Mott insulator (Mott) when the ki- 
netic energy is much smaller than the potential energy and 
a superfluid otherwise. The Mott phase has an excitation gap 
and is incompressible, and therefore, the bosons are localized 
and incoherent, so that a slight change in the chemical po- 
tential does not change the number of bosons on a particu- 
lar lattice site. The superfluid phase, however, is gapless and 
compressible, and the bosons are delocalized and move coher- 
ently. Both of these phases, as well as the transition between 
the two, have been successfully observed with ultracold point- 
like Bose gases loaded into optical lattices M. r9l fiol Hill . 

The on-site BH model takes only the on-site boson-boson 
repulsion into account, i.e. the interaction is short-ranged. A 
more general extended BH model is required when longer- 
ranged interactions are not negligible, e.g. Coulomb or dipole- 
dipole interactions. For instance, an ultracold dipolar Bose 
gas can be realized in many ways with optical lattices Il2ll : 
(ground-state) heteronuclear molecules which have perma- 
nent electric dipole moments, Rydberg atoms which have 
very large induced electric dipole moment, or Chromium-like 
atoms which have large intrinsic magnetic moment, etc. can 
be used to generate sufficiently strong long-ranged dipole- 
dipole interactions. The qualitative phase diagram of this 
model has also been known for a long time II13L IS EE 
E,E>[Hlh an d it has two additional phases: a charge-density 
wave (CDW) as shown in Fig. Q] and a supersolid. Similar to 
the Mott phase, the CDW phase is an insulator with an ex- 



citation gap and it is incompressible. The main difference is 
that an integer number of bosons occupy every lattice site in 
the Mott phase, while the CDW phase has a crystalline order 
in the form of staggered boson numbers (different occupancy 
on different sublattices). As the name suggests, a supersolid 
phase [20], however, has both the superfluid and crystalline 
orders, i. e. both CDW and superfluid phases coexist. There is 
some evidence that this phase exists only in dimensions higher 
than one EE- 



There has been experimental progress in constructing ultra- 
cold dipolar gases of molecules, namely ground-state K-Rb 
molecules, from a mixture of fermionic 40 K and bosonic 87 Rb 
atoms ll2Tll22ll . While this K-Rb is a fermionic molecule, sim- 
ilar principles will allow one to also create bosonic dipolar 
molecules by simply changing the atomic isotopes. Motivated 
by these achievements, in this paper, we analyze the momen- 
tum distribution of the insulating phases of the extended BH 
model, which is the most common probing technique used in 
atomic systems to identify different phases. 



The remainder of this paper is organized as follows. Af- 
ter introducing the model Hamiltonian in Sec. [TH we develop 
two methods in Sec. [ill] to calculate the momentum distribu- 
tion of the insulating (Mott or charge-density-wave) phases 
of the extended Bose-Hubbard model. First we use the ran- 
dom phase approximation (RPA) in Sec. IIII Al and then we 
perform a power series expansion in the hopping t via the 
strong-coupling perturbation theory in Sec. IIII B I The nu- 
merical analysis of the momentum distribution obtained from 
these methods are discussed in Sec. IIII CI and a brief sum- 
mary of our conclusions is presented in Sec. HV] Finally in 
AppendixlAl we comment on some of the issues regarding the 
Wannier functions in the CDW phase. 
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II. EXTENDED BOSE-HUBBARD MODEL 

We consider the following extended BH Hamiltonian with 
on-site and nearest-neighbor boson-boson repulsions 

where iy is the tunneling (or hopping) matrix between sites 
i and j, &j (6^) is the boson creation (annihilation) operator 
at site i, rii — b\bi is the boson number operator, U > is 
the strength of the on-site and Vy is the longer-ranged boson- 
boson repulsion between bosons at sites i and j, and fi is the 
chemical potential. In this paper, we assume is a real sym- 
metric matrix with elements — t for i and j nearest neigh- 
bors and otherwise and similarly for Vij (equal to V > for 
i and j nearest neighbors and zero otherwise), and consider a 
d-dimensional hypercubic lattice with M sites. Note that we 
work on a periodic lattice without an external trapping poten- 
tial. We also assume U > zV where z — 2d is the lattice 
coordination number (number of nearest neighbors). 

The ground-state phase diagram of this model Hamilto- 
nian has been studied extensively in the literature including 
the mean-field fl3tl . quantum Monte Carlo Il4l,[l5ll . density- 
matrix renormalization group flfill . Gutzwiller ansat z IT7i[T8ll . 
and strong-coupling expansion and scaling theory IU9I1 tech- 
niques. When V 0, the ground state now has two types of 
insulating phases. The first one is the Mott phase where, simi- 
lar to the on-site BH model, the ground-state boson occupancy 
is the same for every lattice site, i.e. (rii) = uq. Here, (...) 
is the thermal average, and the average boson occupancy no 
is chosen to minimize the ground-state energy for a given /i. 
The second one is the CDW phase which has crystalline order 
in the form of staggered boson occupancies, i.e. (rii) = n-A 
and (rij) = ns for i and j nearest neighbors. To describe 
the CDW, it is convenient to split the entire lattice into two 
sublattices A and B such that the nearest-neighbor sites be- 
long to a different sublattice. A lattice for which this can be 
done is called a bipartite lattice, and we assume the number 
of lattice sites in each sublattice is the same (M/2). We also 
assume that the boson occupancies of the sublattices A and B 
are ua and ns, respectively, such that ua > The case 
with ua = Ub = no corresponds to the Mott phase. 

When t = 0, it turns out that the chemical potential width of 
all Mott and CDW lobes are U and zV, respectively, and that 
the ground state alternates between the CDW and Mott phases 
as a function of increasing fj, ifll [Tl [Tl [l6i \l% [Til [HI . For 
instance, the ground state is a vacuum (tia = 0,ub = 0) 
for n < 0; it is a CDW with (n A = l,n B = 0) for < 
/i < zV; it is a Mott insulator with (ua = 1, n>B = 1) for 
zV < [i < U + zV; it is a CDW with (n A = 2,n B = 1) 
for U + zV < [i < U + 2zV; it is a Mott insulator with 
(n A = 2,n B = 2) for U + 2zV < [i < 2U + 2zV, and 
so on. As t increases, the range of fi about which the ground 
state is insulating decreases, and the Mott and CDW phases 



dV=0.2U 




dt/U 

FIG. 1: (Color online) We show the chemical potential /i (in units 
of U) versus hopping t (in units of U /d) phase diagram within the 
random-phase approximation for d-dimensional hypercubic lattices 
(it becomes exact for d — * oo). Here the nearest-neighbor repulsion 
scales inversely with d such that dV = 0.217. The red solid line 
correspond to phase boundaries for the Mott insulator to superfluid 
and CDW insulator to supersolid states as obtained from Eq. 



disappear at a critical value of t, beyond which the system 
becomes compressible (superfluid or supersolid) as shown in 

Fig.m 

Identification of these phases in atomic systems loaded into 
optical lattices is a real challenge, and the momentum distri- 
bution of particles has been the most commonly used probing 
technique to distinguish superfluid and Mott phases of the on- 
site BH model. Motivated by these experiments, next we ana- 
lyze the momentum distribution of the insulating phases of the 
extended BH model, paying particular attention to what sig- 
natures one might see that can distinguish the CDW insulating 
phase. 



III. MOMENTUM DISTRIBUTION 

The momentum distribution of the atoms is one of the few 
(and probably the easiest) physical quantity that can be di- 
rectly probed in experiments with ultracold atomic gases. This 
is achieved by time-of-flight absorption imaging of freely ex- 
panding atoms that are released from the trap. Since ultracold 
gases are very dilute, atoms do not interact much with each 
other during this short time-of-flight, and therefore, the parti- 
cle positions in the absorption image are strongly correlated 
with their velocity distribution given by their momentum dis- 
tribution at the moment of release from the trap. 

The momentum distribution n(k) is also easy to calculate, 
and it is defined as the Fourier transform of the one-particle 
density matrix p(r, r') = (^^(r)?/'(r')), such that 

n(k) = [ dr [ dY' P {r,v'y^ ir - r '\ (2) 



3 



where ip >(r) [tp(r)] is the boson creation (annihilation) field 
operator, and k is the momentum. We expand the field oper- 
ators in the basis set of Wannier functions such that ip{r) = 
(1/a/M) W(r — He)bt, where M is the number of lattice 
sites, and the Wannier function W(r — R^) is localized at site 
£ with position R^. Here the summation index £ 6 {A, B} 
includes the entire lattice. 

In this paper, we use two methods to calculate the mo- 
mentum distribution of the insulating phases of the extended 
BH model. First we calculate n(k) via the RPA theory in 
Sec. IIII Al and its result corresponds to the exact result for 
the infinite-dimensional limit. Then, in Sec. HUB I we calcu- 
late n(k) as a power series expansion in the hopping t via 
the strong-coupling perturbation theory. We also verify that 
our strong-coupling expansion recovers the RPA result in the 
infinite-dimensional limit when the latter is expanded out in t 
to the same order. This provides an independent cross-check 
of the algebra as discussed next in detail. 

A. Random Phase Approximation (RPA) 

Using the standard-basis operator method developed by Ha- 
ley and Erdos B23I1 . and following the recent works on the 
on-site BH model JM IH IH, Hi 0], here we obtain the 
equation of motion for the insulating phases of the extended 
BH model. This approximation is a well-defined linear opera- 
tion in which thermal averages of products of operators are re- 
placed by the product of their thermal averages. In accordance 
with this approximation, the three-operator Green's functions 
are reduced to two-operator ones B23I1 . Therefore, the RPA 
method allows us to calculate the single-particle Green's func- 
tion G(k, iuj n ) = —(i{)(\s.,iu; n )ipi(k.,i(j n )) in momentum (k) 
and Matsubara frequency (iuj n ) space, from which the spec- 
tral function A(k,ui) = — (l/7r)ImG(k, iu) n — > uj+ie) canbe 
extracted by analytical continuation. Here the angular brack- 
ets denote the standard trace over the density matrix. Notice 
that the spectral function should always satisfy the sum rule 
A(k, cu)duj = 1, due to the bosonic commutation rela- 
tions of the creation and annihilation operators in the Heisen- 
berg picture at equal times. Then the momentum distribution 
n(k) = (ip< (k)^(k)) (at zero temperature) can be easily ob- 
tained from the spectral function n(k) = — f_ A(k, u))cku, 
i.e. 

1 f° 

n(k) = — ImG(k, iuj n — ► u> + ie)du, (3) 

1" J-oo 

which measures the spectral weight of the hole excitation 
spectrum. 

Expanding the field operators given in Eq. (f2| in the basis 
set of Wannier functions, the momentum distribution becomes 

„(k) = i^M^^^^e-*-^-^), (4) 

where W(k) = J driy(r)e lk r is the Fourier transform of 
W(r). Here the summation indices I e {A, B} and t' G 



{A, B} include the entire lattice. Since VK(k) is a nonuniver- 
sal property of the lattice potential, and it has nothing to do 
with the extended BH model on a discrete periodic lattice, we 
ignore this function in this paper by setting it to unity. 

But before beginning the discussion of our formal treatment 
of the theory, we want to comment further on the subtle fea- 
tures that arise for the momentum distribution n(k) in an or- 
dered phase, when the lattice periodicity is further broken by 
the spontaneous appearance of the CDW phase with a lower 
lattice periodicity. This system becomes that of a lattice with 
a basis, as the A and B sublattices now have a different oc- 
cupancies of particles on them. When examining n(k) on the 
lattice, we evaluate the one-particle density matrix at each lat- 
tice site p(r,r') — * /9(R,,Rj) = (see the appendix for 
a further discussion of how one goes from the continuum to 
the lattice and how Wannier functions enter into the calcula- 
tion). The integral in Eq. © is replaced by a summation that 
extends over all lattice sites of the original lattice (before the 
CDW order occurred). We can break this summation up into 
terms that involve solely the A sublattice, solely the B sub- 
lattice, and terms that mix the A and B sublattices. One can 
immediately see that the terms restricted to one of the sublat- 
tices are periodic with the periodicity of the reduced Brillouin 
zone, while the mixed terms are only periodic with respect to 
the full Brillouin zone. If we assume the Wannier functions 
are identical for the A and B sublattices, then this uniform 
weighting of the different contributions yields the correct mo- 
mentum distribution; in general, one potentially has different 
weightings of the three different components. A full discus- 
sion of this issue is beyond this work, where we focus on the 
properties of the pure discrete lattice system, not on the exper- 
imental systems which have the additional real-space structure 
arising from the spatial continuum. 

The fluctuations are not fully taken into account in the RPA 
method, however it goes beyond the mean-field approxima- 
tion for low-dimensional systems, and it becomes exact for 
infinite-dimensional bosonic systems recovering the mean- 
field theory. The RPA method has recently been applied to 
describe the superfluid and Mott phases of the on-site BH 
model rf25l l26ll . and its results showed good agreement with 
the experiments. Motivated by these earlier works, here we 
generalize this method to describe the insulating phases of the 
extended BH model. 

Keeping in mind our two-sublattice system, the single- 
particle Green's function in momentum and frequency space 
can be written as G(k, iu) n ) = (1/2) J2s S' Gss'fe, iu> n ), 
where the indices S and S' label sublattices {A, B}, and 
G ss ,(k,Mn) = (yM)^ ieS t, eS ,G M ,(iu„)e- ik < R *- R ^ 
is the Fourier transform. Here the summation indices I 6 A or 
B and i' E A or B include only one sublattice and the Green's 
function is defined only at the different lattice positions. Since 
there are M/2 lattice sites in one sublattice, a factor of 2 
appears in this expression. Note that this is just a rewriting 
of the summation over all lattice sites that explicitly shows 
the contributions from the different sublattices. The RPA 
equations have the following form in position and frequency 
space Gw(iuj n ) = G°(iu n ) \b w + J2e>> Ju"Gi"e'(wn)] , 
where G°(icj„) and J«» are given below Eq. ©. Here 
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the indices £, I 1 and I" e {^4,5} include the entire lat- 
tice. Using the Fourier transforms, the RPA equation in 
momentum and frequency space becomes Gss' (k, iuj n ) — 
G a s (k,ioj n ) [5 S s> + J2s" Jss>> (k)G s >> s>(k, iu> n )} ■ This ex- 



pression defines a set of coupled equations for the 
functions GaaQl, iui n ), G AB (k, iu> n ), G B A(k,iw n ) and 
G BB (k, iu> n ). These equations can be easily solved to obtain 



G(k,iuj n ) = 



[G° A (k,iw n ) + G B (k,uu n )]/2 + [Jab (k) + J BA (k) - JaaM ~ J B B(k)]G A (k, iw n )G B (k, kj n )/2 



1 - J A A{k)G A {k,iu n ) - J BB (k)G B (k,; 



[Jab (k) Jjja (k) - J AA {k) J BB (k)]G A (k, iu n )G° B (k, 



VjJ n 



which corresponds to the general single-particle Green's function within the RPA. 



(5) 



We use Eq. (0 to obtain the single-particle Green's function 
for the insulating (Mott and CDW) phases of the extended BH 
model. Since hopping is allowed between nearest-neighbor 
sites that belong to different sublattices, in Eq. (O we find 
JaaW = J BB (k) = and J AB (k) = J BA (k) = e(k), 
where e(k) is the Fourier transform of the hopping matrix 
(also called the band structure). For e?-dimensional hyper- 
cubic lattices considered in this paper, the energy dispersion 
becomes e(k) = — 2t ^ i=1 cos(fcia), where a is the lattice 
spacing. This then yields the following expression for the 
Green's function: 



Gi ns (k,iu n ) 



GjJMd + e(k)G A (im n )G B (iLo n ) 
l-£ 2 (k)G»K)G°K) ' 



(6) 



tion leads to 



H 2 



( n A + 1 ^ ua_ 



n B 



ppar 
B 



n B 

B 



(9) 



which is a quartic equation for /i, and it coincides with our 
earlier result lfl9ll . Notice that Eq. (O reduces to the usual 
expression for the phase boundary of the on-site BH model 
when riA — n B = no and V = 0. Having discussed the gen- 
eral RPA formalism for the insulating phases of the extended 
BH model, next we analyze the momentum distribution of the 
Mott and CDW phases separately. 



1. Mott Phase 



where G° vr (^n) = [G A {iu n )+G B {iu n )]/2. The k- 
independent functions G A (ioj n ) and G B (iui n ) correspond to 
the single-particle local Green's functions for sublattices A 
and B, respectively, at zeroth order in t. They have the famil- 
iar form 



G° A (ico n ) = 
G° B (iuj n ) = 



ua + 1 



riA 



j-ipar 
^A 



^hol : 



n B + 1 



n B 



^hol : 



(7) 
(8) 



where E A 



Un A + zVn B — ll and E B m — Un B + zVriA — 
ix are the zeroth-order particle excitation spectrum in t (the 
energy required to add one extra particle) for sublattices A and 
B, respectively, and similarly E A ° l — ~U(riA — l) — zVn B + 
li and E B ° l = — U(n B — 1) — zVtia + Li are the zeroth-order 
hole excitation spectrum in t (the energy required to remove 
one particle). Notice that Gi ns (k, iw n ) = G° vr (iw n ) at zeroth 
order in t, as one may expect. 

The poles of Gi ns (k, iuj n ), the condition 1 = 

e 2 (k)G A (iuj n )G B (iuj n ), give the k-dependence of the par- 
ticle and hole excitation spectrum. The insulating phase be- 
comes unstable against superfluidity when any of the excita- 
tion energies becomes negative at k = 0. In addition, the 
poles of Gi ns (k,iuj n ) at (k = 0,iu n — 0), i.e. the con- 
dition 1 = e 2 (O)G^(0)G s (0), gives the mean-field phase 
boundary between the incompressible (Mott or CDW) and the 
compressible (superfluid or supersolid) phases. This condi- 



The single-particle Green's function for the Mott phase can 
be obtained from Eq. (O by setting n A — n B — n . This 
leads to 



Gmou (k, iio n 



G^iuJn) 



1 - e(k)G° (iu> v 



(10) 



which has the same form with that of the Green's function 
of the Mott phase in the on-site BH model JH [H, |28ll . 
Here, G A (iuj n ) = G B (iuj n ) = Go(£cJ ra ). The function 



GMott(k i^r, 

'(k), 



-E^i 



has two poles at iw n = Eq 3,1 (k) and ioj n 



E^(k) = E*™ -[U- s(k) - E (k)] /2, 



4 lol (k) 



E 



hoi 



[U + e(k) - E (k)] /2, 



(11) 
(12) 



corresponding to the particle (the energy required to add 
one extra particle) and hole (the energy required to re- 
move one particle) excitation spectrum, respectively, where 
E (k) = v / e 2 (k) + 2U(2n + l)e(k) + U 2 . Notice that the 
Mott insulator becomes unstable against superfluidity when 
E$ ar (0) = or E W \0) = °. and these conditions co- 
incide with the mean-field condition given in Eq. (0 when 
n A =n B = n . 

Therefore, the Green's function for the Mott phase can be 
written as 



GMott(k i^n 



G par (k) 



Co ( k ) 



E^(k) 



E$°\k) 



, (13) 
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where the coefficients (or the spectral weights) are functions 
of the excitation spectrum 



'00 = h ° 



C hol (k) = 



r (k) + Uno 



Eft* 



E pa 



+ E$°\k) 
Un 



Eq 



E pi 
^0 



r (k)+E$°\k) 



(14) 
(15) 



perturbation theory up to second order in t. Notice that the 
momentum distribution is flat and equals the average filling 
fraction riMott (k) = no at zeroth order in t, corresponding to 
vanishing site-to-site correlations. 



CDW Phase 



Using the definition given above Eq. (0, the spectral function 
for the Mott phase can be easily obtained from Eq. ( fT3b . lead- 
ingto^Mott(k,^) = C pm \k)S[uj-E pai (k)}+C^ 1 (k)S[tu+ 
Eq (k)], where 5 (x) is the Delta function defined by S(x) = 
(1/ir) lim e ^o e / (x 2 + e2 )- Notice that this function satisfies 
the sum rule mentioned above Eq. (0, since the coefficients 
satisfy Cg al (k) + C^ 1 (k) = 1. The momentum distribution 
measures the spectral weight of the hole excitation spectrum 
as defined in Eq. (0, and for the Mott phase it is given by 



"Mott(k) = -C hol (k) 



U(2n + 1) + e(k) _ I 
2E (k) r 



(16) 



which is identical to the riMott(k) of the on-site BH 
model 1251, l26ll . Therefore, at the RPA level, UMott(k) is in- 
dependent of V which is mainly because of the underlying 
mean-field Hamiltonian that is used in the RPA formalism (we 
remind that fluctuations are not fully taken into account within 
RPA). For instance, the mean-field phase boundary condition 
given in Eq. (0 shows that the Mott lobes are separated by 
zV, but their shapes and, in particular, the critical points are 
independent of V. This point will become more clear in 
Sec. IIIIBI where we analyze n(k) via the strong-coupling 



In contrast to the Green's function of the Mott phase, 
the single-particle Green's function for the CDW phase 
GcDw(k, ioj n ) has four poles. Two of them correspond to 
the particle and the other two to the hole excitation spectrum 
of sublattices A and B. Unfortunately, general expressions 
for these poles are not analytically tractable since the condi- 
tion 1 = e 2 (k)G A (iuj n )G B (iuj n ) defines a quartic equation 
for iu> n ; they can be easily obtained numerically for any given 
CDW lobe as shown in Sec. IIIICI Assuming that the exci- 
tation spectrum is known, the Green's function for the CDW 
phase can be written as 



G 



CDW 



(k, iWr, 



^par 



(k) 



C\°\k) 



E p /\k) 



(k) 



+ 



(k) 



(k) 



£T(k) 



E%°\ky 



(17) 



where _E^ ar (k) and E B m (k) are the particle (the energy re- 
quired to add one extra particle) and E 1 ™ 1 (k) and Eg o1 (k) are 
the hole (the energy required to remove one particle) excita- 
tion spectrum. The coefficients (or the spectral weights) are 
functions of the excitation spectrum, such that 



^par 



(k) = 
'(k) = 



^ A 



(k) 



Choi 
B 



A>(k) 



[E*?(k) 
D (k) 



^D 1 (k)E^(k) 
■E^(k)][E^(k) 



D 2 {k)[ET(k)f 



- Dl (k)E^(k) 



[E^(k) 
D (k) 



-E p /\k)][E 



[^ ar (k)] 3 

^(k)]K ar (k)+^ ol (k)]' 
D 2 (k)[E B ^(k)r + [E^(k) 



<(k)+E A °\k)][E B ar (k) + E B °\k)} 



D 1 (k)EY(k) + D 2 (k)[E^ (k)] 2 - [E A 



B 

,hol (k)] 3 



(k) = 



D (k) 



- i^oi (k) phoi (k) + E ^(k)]{E\°\k) + E p /*(k)} ' 
D^E^Hk) 



D 2 (k){E B °\k)f 



I^ ol (k)] 3 



[E\°\k) - B^(k)][^ ol (k) + ^r(k)][^ ol (k) + E p / r (k 



(18) 
(19) 
(20) 
(21) 



Here, the coefficients £>o(k), Di(k), and D2(k) are func- 
tions of the zeroth-order excitation spectrum in t defined be- 
low Eq. (0, and are given by 



A>(k) = - 
(Un A 



hoi 



[ETE A 



(Un B + E 



hoh 



ETE B 



hoi 



E A ° 1 )] /2 + e(k)(Un A + E A ol )(Un B 



(22) 



£>i(k) = [(E A ° l - E p / T )(Un B + E B ° l ) + (E™ - E^) 
{Un A + EY) - E p / r E A ° l - E p ™E B ° l ] /2 
+ e(k)(Un A + Un B + E A ° l +E B ° l ), (23) 

and 

D 2 (k) = (Un A + Un B - E paI - E pal ) /2 + E A ° l + E B ° l 
+ e(k). (24) 

Using the definition given above Eq. <[3j, the spec- 
tral function for the CDW phase can be easily obtained 
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from Eq. ( [T7l i, leading to j4cDw(k, u) = C^ ar (k)S[uj — 
E^(k)}+Cr(m^+E\ \k)]+C p -(k)S[co-E^(k)] + 
C B ol (k)S[cj + £^ ol (k)]. Notice that this function satisfies the 
sum rule mentioned above Eq. (O, since the coefficients sat- 
isfy C^ ar (k) + C^ ol (k) + C p B a \k) + C B ol (k) = 1. The mo- 
mentum distribution measures the spectral weight of the hole 
excitation spectrum as defined in Eq. ([3j, and for the CDW 
phase it is given by 



ncDw(k) = ~C\°\k)-C h B °\k). 



(25) 



This expression has a highly non-trivial dependence on t, and 
it has to be solved numerically together with the excitation 
spectrum. However, it can be analytically shown that the mo- 
mentum distribution is flat and equals the average filling frac- 
tion ncDw(k) = (n A + ns)/2 at zeroth order in t, corre- 
sponding to vanishing site-to-site correlations. To provide an 
independent check of the algebra (and to extend to finite di- 
mensions), we next calculate n(k) as a power series expan- 
sion in the hopping t via the exact strong-coupling perturba- 
tion theory in d dimensions. 



B. Strong-coupling Perturbation Theory 

To determine the momentum distribution of the insulat- 
ing phases, we need the wavefunction of the insulating state 
l^ins) as a function of t. We use the many-body version of 
Rayleigh-Schrodinger perturbation theory in the kinetic en- 
ergy term ll29ll to perform the expansion (in powers of t) for 
l^ins) needed to carry out our analysis. A similar expansion 
for the ground-state energies was previously used to discuss 
the phase diagram of the on-site BH model iH [4 ], and it has 
recently been applied to the extended BH model 111911 . For the 
on-site BH model, extrapolated results of these expansions 
showed an excellent agreement with recent quantum Monte 
Carlo simulations JH |g]. A high-order strong-coupling ex- 
pansion for the ground-state energies has now been extended 
to all dimensions and fillings M30H . and a high-order expansion 
for the wavefunction has also been used to describe the Mott 
phase in one-dimensional systems l3~ill . 

For our purpose, we first need the ground-state wavefun- 
tions of the Mott and CDW phases when t = 0. To zeroth 
order in t, the insulator (Mott or CDW) wavefunction can be 
written as 



(OK 



M/2 

n 




(26) 



state (here, we remind that the lattice is divided equally into 
A and B sublattices). In principle, we can apply the pertur- 
bation theory on l^w) to calculate l^ins) up to the desired 
order. However, since the number of intermediate states in- 
creases dramatically due to the presence of nearest-neighbor 
interactions, we perform this expansion only up to second or- 
der in t. The (unnormalized) wavefunction for the insulating 
state can then be written as 



Ins 



I* 



(Oh 



Ins/ 



+ 



E 



E 

T / 



tt" I Ins/ 



— \<l)+0(t 3 ), (27) 



{m',m}^|*'°>> 



Eom> Eq 



where T m0 = - J2 S ,S' T,ees,£'eS' *«'( m l b ]M*W is the 
hopping matrix element between the first-order intermediate 

state \m) and the zeroth-order state |^jns)> an d T mm i is be- 
tween |m) and the second-order intermediate state \m'), and 

E 0m — E^l — Em . Here the summation indices I S {^4, B} 
and £' G {A, B} include the entire lattice, and S and S' label 
sublattices {A,B}. The \m) states are connected to l*^) 
state with a single hopping, and similarly \m') states are con- 
nected to | to) states with a single hopping. However, the \m') 
state must be different from the l*^) state. 

To calculate the momentum distribution, we need the 
normalized wavefunction for the insulating state |Wi ns ) = 
l^ins)/ yj (V'inslV'ins), where the normalization up to second 
order in t is given by 



(■0Im#Iii 



n A (n B + l)Mzt 2 /2 



[U (tia - n>B - 1) + V(zn B - zn A + 1)Y 
n B (riA + l)Mzt 2 /2 



[U(riB - tia - 1) + V(zn A - zjib + 1)]" 



o(t 4 ). 



(28) 



Here, 

(ml* 



= 2d is the lattice coordination number. Since 



(Oh 



(m"m 



(oh 



0, the first and third 



where M is the number of lattice sites, and 10) is the vacuum 



order terms in t vanish in the normalization. In general, all 
odd-order terms in t vanish. 

A lengthy but straightforward calculation leads to the 
momentum distribution, [defined in Eq. n(k) = 

(l/M)J2 ee ,(y lns \blb e ,\y lns )e~* k < R '- R t'\ up to second 
order in t as 
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nias(k) 



n A + n B 



U A (n B + 1) 



n B {n A + 1) 



U(n,A — n B — 1) + V(zn B - zn A + 1) U(n B — n A - 1) + V(zn A — zn B + 1) 
n A (n B + l) n B (n A + l) 



e(k) 



■ 2 [/7(n A - nu - 1) + y(zn s - zii A + l)] 2 2 [L/(n B - n A - 1) + ^(ot^ - zn B + l)] 2 



+ 1) 



n B (n A + 1) 



[7 [Z7(>M - n B - 1) + V(zn B - zn A + 1)] ?7 [U (n B - n A - I) + V(zn A - zn B + 1)] 
(n A +n B + l) [e 2 (k) - 2dt 2 ] + 0(t 3 ). 



(29) 



In the definition of the momentum distribution, the summation 
indices £ S {A, B} and £' e {A, _B} include the entire lat- 
tice. Here e(k) = -(2/M) Etes^'eS' i«'e lk - (Rf - Rf,) 
is the Fourier transform of the hopping matrix 
t^/ (energy dispersion), and £ 2 (k) — 2dt 2 = 

(2/M)Y /{e j ll}eS j'eS' t u' t e'e" e ' kiRe " Re "^ where the 
summation indices {£, £"} e A (or B) and £' € B (or A) 
include only one sublattice. Since there are M/2 lattice sites 
in one sublattice, a factor of 2 appears in these expressions. 
To zeroth order in t, Eq. (|29l shows that ni ns (k) is flat and 
equals the average filling fraction (n A + n B )/2. However, it 
develops a peak around k = and a minimum around k = n 
at first order in t. These general observations are consistent 
with the RPA results shown in Eqs. ([Tol l and ( [25l >. 

Equation ( |29l is valid for the insulating phases of all d- 
dimensional hypercubic lattices. For instance, when n A — 
n>B = no, Eq. (l29l reduces to the momentum distribution for 
the Mott phase, i.e. 

e(k) 

n-Mott(k) = n - 2n (n + 1) + n (n + 1) 

q/t _ 2V 

(2n + 1) [ £ 2 (k) - 2dt 2 } u{u _ v)2 + 0(t 3 ). (30) 

This expression recovers the known result for the on-site BH 
model when V = IT321 l33ll . In addition, in the d — > oo limit, 
we checked that Eqs. ( 1291 and ( 1301 agree with the RPA solu- 
tions (which are exact in this limit) given in Eqs. d25l l and ( [Tol 
when the latter are expanded out to second order in t, provid- 
ing an independent check of the algebra. One must note that 
the terms 2V and V that appear in the numerator and denom- 
inator of Eq. d30l vanish in the limit when d — * oo because 
V oc 1/d. Next we compare the RPA results with those of the 
strong-coupling perturbation theory. 



C. Numerical Results 

Since the momentum distribution of the CDW phase given 
in Eq. ( fZSb has a highly nontrivial dependence on t, it has to 
be solved numerically together with the excitation spectrum. 
Next we set dV — 0.2U and solve this equation for the first 
CDW lobe. For this parameter, we remind that the t = 
chemical potential width of all Mott and CDW lobes are U 
and OAU, respectively, and that the ground state alternates 



between the CDW and Mott phases as a function of /i. For 
instance, the ground state is a vacuum (no = 0) for fi < 0; it 
is a CDW with (n A = 1, n B = 0) for < fi < OAU; it is 
a Mott insulator with (n = 1) for OAU < \i < 1AU; it is a 
CDW with (n A = 2, n B = 1) for 1.417 < \i < 1.8U; it is a 
Mott insulator with (n = 2) for 1.8U <H< 2.8U. 

In Fig. |2] the results of the RPA calculation given in Eq. d25l ) 
are compared to those of the second-order strong-coupling 
perturbation theory given in Eq. d29l for a (d — 2)- and 
(d — > oo) -dimensional hypercubic lattices. In this figure, we 
show the momentum distribution «cdw (k) as a function of 
e(k) / (dt) for two sets of parameters. In Fig. |2|a), we choose 
dt = 0.05U and \i = 0.2U which approximately corresponds 
to the center of the first CDW lobe. For this parameter set, 
deep inside the CDW lobe, the momentum distribution has 
a peak at e(k) = — 2dt corresponding to the k = point, 
and it has a minimum at e(k) = 2dt corresponding to the 
k = (w, 7r, . . •) point. This is very similar to what hap- 
pens in the Mott phase. However, in Fig. |2fb), we choose 
dt sa 0.083J7 and n w 0.16C/ which approximately corre- 
sponds to the tip of the first CDW lobe. For this parameter set, 
close to the CDW-supersolid phase transition, the momentum 
distribution has two peaks: a large peak at e(k) = —2dt corre- 
sponding to the k = point, and a smaller one at e(k) = 2dt 
corresponding to the k = (%,%,...) point. The second peak 
is unique to the CDW phase and it does not occur in a Mott 
phase. Notice that both the RPA and second-order strong- 
coupling expansion give qualitatively similar results (although 
the peak is much sharper and has lower weight in the exact so- 
lution). 

One might have expected to always see the peak in the mo- 
mentum distribution at the k = (tt, tt, . . .) point due to the 
reduced periodicity of the CDW order. But because the mo- 
mentum distribution involves four terms corresponding the the 
AA, AB, BA, and BB sublattice combinations, only the first 
and last terms are periodic in the reduced Brillouin zone (see 
our discussion given in the appendix). Deep inside the CDW 
lobe, the presence of a large gap in the one-particle excitation 
spectrum produces an exponential decay of the one-particle 
correlations which suppresses this peak in the momentum dis- 
tribution as can be seen in Fig. |2ja) (this point has already 
been discussed in Ref. 1341) . This essentially occurs because 
there is a cancellation of the peak that arises from the AA 
and BB contributions with the results from the AB and BA 
pieces, similar to what happens in the Mott phase. However, 
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(a) dt=0.05U and u=dV=0.2U 




0.2 1 1 1 1 1 

-2-1012 

8(K)/(dt) 

(b) dt=0.083U, dV=0.2U and u=0.16U 




-2-1012 

e(k)/(dt) 



FIG. 2: (Color online) Momentum distribution ncDw(k) versus 
e(k)/(dt) for a (d = 2)- and (d — > oo)-dimensional hypercubic lat- 
tices. Panel (a) has the nearest-neighbor boson-repulsion satisfying 
dV = 0.2(7, the hopping satisfying dt — 0.05(7 and the chemical 
potential set by fi = 0.2(7 corresponding approximately to the center 
of the first CDW lobe, and panel (b) has dV = 0.2(7, dt ps 0.083C/, 
and fi ~ 0.16 corresponding approximately to the tip of the first 
CDW lobe. The solid (red) lines correspond to the RPA, and the 
dashed and circled lines to the second-order strong-coupling pertur- 
bation theory for different dimensions. The peak occurs at the zone 
corner only when the hopping is close to the tip of the CDW lobe. 



close to the tip of the CDW lobe, the peak emerges in the exact 
solution of the RPA as shown in Fig.|2|b). To some extent, this 
peak also emerges in the solutions of the second-order strong- 
coupling perturbation theory. Notice that the peak is under- 
emphasized in the strong-coupling theory since the theory is 
exact only deep inside the CDW lobe, and it becomes quanti- 
tatively inaccurate for large values of dt/U close to the tip of 
the CDW lobe. We remark that an unphysical peak appears at 
k = (tv, 7T, . . .) in the strong-coupling perturbation theory for 
the Mott phase (not shown), which signals the breakdown of 
the second-order expansion. 

As a further check of the accuracy of our second-order 



strong-coupling expansion, in Fig. [3] we compare the d = 2 
and d — > oo limits of Eq. d29l ) to the RPA method given 
in Eq. ( f25T > which corresponds to the exact solution in the 
latter limit. In this figure, we show ncD\v(k = 0) an d 
^CDw(k = 7r) as a function of dt/U when fi = 0.2(7. In 
d — 2 dimensions, the RPA and second-order strong-coupling 
expansion gives qualitatively similar results for small values 
of dt/U, i.e. deep inside the CDW lobe. However, in the 
d — ► oo limit, the results of the RPA and the second-order 
strong-coupling expansion match exactly for small values of 
dt/U (as they must). Close to the tip of the CDW lobe, the 
RPA and strong-coupling results differ substantially from each 
other signalling the breakdown of the second-order expansion. 
However, both theories show that ticdw(O) is an increasing 
function of dt/U as one may expect. This is because the range 
of /i about which the ground state is a CDW decreases as dt/U 
increases from zero, and the CDW phase become a supersolid 
at a critical value of dt c ~ Q.08U. Beyond this point, n(0) di- 
verges due to the appearance of a condensate, corresponding 
to the macroscopic occupation of the k = state. 

Note that we do not attempt to perform a scaling analysis 
of the momentum distribution for the CDW phase. The rea- 
sons why are twofold. First, we only have the series through 
second order, which probably is too short to be able to prop- 
erly fit to a phenomenological scaling form, and second, we 
cannot extract the analytic scaling form from the RPA calcu- 
lation anymore, so guessing an appropriate phenomenological 
form has less guidance than for the Mott phase. A scaled the- 
ory would be expected to be accurate for all values of t within 
the insulating phases, as has been recently shown for the Mott 
phase of the on-site BH model 11320 . 



H=dv=0.2U 




0.02 0.04 0.06 0.08 
dt/U 



FIG. 3: (Color online) Momentum ditributions at specific momen- 
tum points 7icDw(k = 0) and ncDw(k = n) versus dt/U for 
(d — 2)- and (d — ► oo) -dimensional hypercubic lattices. The chem- 
ical potential /i = dV corresponds to the first CDW lobe, and the 
nearest-neighbor repulsion is set to dV = 0.2(7. The solid line cor- 
responds to the RPA and the dashed and circled lines to the second- 
order strong-coupling perturbation theory for different dimensions. 



9 



IV. CONCLUSIONS 

We developed two methods to calculate the momentum 
distribution of the insulating (Mott and charge-density-wave) 
phases of the extended Bose-Hubbard model with on-site and 
nearest-neighbor boson-boson repulsions on d-dimensional 
hypercubic lattices. First we analyzed the momentum distri- 
bution within the random phase approximation, which cor- 
responds to the exact solution for the infinite-dimensional 
limit. Then we used the many-body version of the Rayleigh- 
Schrodinger perturbation theory in the kinetic-energy term, 
and derived the wavefunction for the insulating phases as a 
power series in the hopping t, to calculate the momentum dis- 
tribution via the strong-coupling perturbation theory. A sim- 
ilar strong-coupling expansion for the ground-state energies 
was previously used to discuss the phase diagram of the on- 
site BH model |Hl 01, and it has recently been applied to the 
extended BH model fl9ll . 

The agreement between the second-order strong-coupling 
expansion and that of RPA method is only qualitative in low- 
dimensional systems. This is not surprising since the fluctu- 
ations are not fully taken into account in the RPA method. 
However, we showed that our strong-coupling expansion 
matches exactly the RPA result (as it must) in the infinite- 
dimensional limit when the latter is expanded out in t to the 
same order. We believe some of these results could poten- 
tially be tested with ultracold dipolar Bose gases loaded into 
optical lattices. This work can be extended in several ways 
if desired. For instance, one could calculate the momentum 
distribution up to third order in i, and develop a scaling theory 
with the help of the RPA results (or a good phenomenologi- 
cal guess for the scaling form of the momentum distribution). 
The scaled theory is expected to be accurate for all values of 
t within the insulating phases, as has been recently shown for 
the Mott phase of the on-site BH model ll32tl . 
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APPENDIX A: EFFECTIVE CDW HAMILTONIAN 

In this Appendix, we comment on some of the subtle is- 
sues regarding Wannier functions in the CDW phase. When 
the particle occupancies show a CDW order, we can think of 
the combination of the CDW order plus the lattice potential as 
an effective lattice potential such that the effective potential 
is different for each sublattice. In other words, CDW order 
creates an effective potential which depends on the particle 
occupation of the sublattice. In fact, having different effective 
lattice potentials on two sublattices could be thought of as the 
reason for having a CDW order at the first place. Equivalently, 
this is like considering the mean-field Hamiltonian with CDW 



order as the starting point for determining the Wannier wave- 
functions, with the symmetry explicitly broken between the A 
and B sublattices. 

This observation suggests that in contrast to the Mott phase 
where all lattice sites are identical and the Wannier func- 
tions are exactly the same for both sublattices, i.e. Wa{y) = 
Wb(y) — Wo(r), the Wannier functions depend on the sub- 
lattice when the CDW order exists, i.e. Wa{f) ^ Ws(r). 
Throughout this paper, we assume that the Wannier func- 
tions are equal (or at least similar) in sublattices A and B. 
However, depending on the CDW order (e.g. tia S> ns) 
and the lattice potential, the Wannier functions of one sub- 
lattice may become substantially different from that of the 
other. In such a case, the field operator can be expanded 

as V(r) = (1/VM) J2s E/ 6 s Wsi ? ~ R <) 6 <> where M is 
the number of lattice sites, S labels sublattices {A, B}, and 
M / 5(k) = J drIVs(r)e lk ' r is the Fourier transform. Here the 
summation index I S {A, B} includes the entire lattice. 

When WOt(r) ^ Ws(r), the strength of the on-site boson- 
boson repulsion also depends on the sublattice, since the ef- 
fective interaction Ug S = g J dr\ Ws(r)| 4 is larger for deeper 
potentials, where g is the bare boson-boson repulsion of the 
continuum Hamiltonian. Therefore, the effective Hamiltonian 
that describes the CDW phase can be written as 

^Dw = ^ cff E 5 ^-^ E ^ 

ie{A,B} 

+ -A. rii {% - 1) + -f- E % (% - !) 

+ Vab E (A1) 

(ieA,jeB) 

where the notation (i,j) corresponds to nearest-neighbors. 
Here the effective hopping element t A s B — t c g A = t be- 
tween the two sublattices is given by t = — J drW' A (r — 
R-i) [- V 2 /(2m) + Vol (r)]Ws (r - Rj), where m is the mass 
of particles, and Vol (r) is the lattice potential, and V A % = 
g J dr\WA{r - Ri)| 2 |W B (r - Rj)| 2 is the effective nearest- 
neighbor boson-boson repulsion. 

When W^4(r) ^ W^(r), the momentum distribution given 
in Eq. (01 becomes 

n CD w(k) = ^E E W* s (k)W s >(k) 

s,s' ees,e'eS' 

(6|V)e- lk -( Rf - R ^, (A2) 

since the boson creation and annihilation operators have dif- 
ferent weights depending on their acting sublattice. Here the 
summation indices £ E {A, B} and £' 6 {A, B} include 
the entire lattice. This summation breaks up into terms that 
involve solely the A sublattice, solely the B sublattice, and 
terms that mix the A and B sublattices. One can immediately 
see that the terms restricted to one of the sublattices are peri- 
odic with the periodicity of the reduced Brillouin zone, while 
the mixed terms are only periodic with respect to the full Bril- 
louin zone. In general, these terms have different weightings 
when Wannier functions differ on two sublattices. A detailed 
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analysis of the CDW Hamiltonian given in Eq. (IAU and its they will be addressed elsewhere, 
momentum distribution is beyond the scope of this paper and 
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